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Abstract 

We present here a new image inpainting algorithm based on the Sobolev 
gradient method in conjunction with the Navier-Stokes model. The orig- 
inal model of Bertalmi'o et al. [Proc. Conf. Comp. Vision Pattern Rec, 
2001] is reformulated as a variational principle based on the minimization 
of a well chosen functional by a steepest descent method using Sobolev 
gradients. This new theoretical framework offers an alternative of the di- 
rect solving of a high-order partial differential equation, with the practical 
advantage of an easier and more flexible computer implementation. In 
particular, the proposed algorithm does not require any constant tuning, 
nor advanced knowledge of numerical methods for Navier-Stokes equa- 
tions (slope limiters, dynamic relaxation for Poisson equation, anisotropic 
diffusion steps, etc). Using a straightforward finite difference implemen- 
tation, we demonstrate, through various examples for image inpainting 
and image interpolation, that the novel algorithm is faster than the orig- 
inal implementation of the Navier-Stokes model, while providing results 
of similar quality. 

The paper also provides the mathematical theory for the analysis of the 
algorithm. Using an evolution equation in an infinite dimensional setting 
we obtain global existence and uniqueness results as well as the existence 
of an w-limit. This formalism is of more general interest and could be ap- 
plied to other image processing models based on variational formulations. 



1 Introduction 

Image inpainting refers to the process of filling in an occluded region in an image 
so that the seam between the original image and the inpainted region is unde- 
tectable by a typical viewer. Some examples of applications are restoring old 
photographs, removing unwanted objects from images such as text, or removing 
certain features of an image such as eye glasses. Image inpainting has been 
traditionally the work of professional artists. However, in recent years, much 
work has been done in the area of digital image inpainting. The idea of digital 
image inpainting is that the algorithm, using information from elsewhere in the 
digital image, fills in the unknown region so that a typical viewer cannot judge 
what parts of the image are from the original image and what parts are painted 
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in. In general, the area to be altered consists of many pixels, complex geome- 
tries, and possibly spans the entire image. This establishes image inpainting as 
a challenging computational problem. 

Currently there are a number of models and algorithms that are used in 
digital inpainting. The pioneering works in digital image inpainting include [1] , 
[TT] , [in] , and [3D] . In [3D] , the authors present a method where the region to be 
inpainted is viewed as an occluded area. Then one connects T-junctions with 
the same grayscale values at the boundary using elastica minimizing curves. In 
[T5] . the authors extend on this idea by using a variational formulation. They 
inpaint the region with unknown pixel values by connecting isophotes at the 
boundary using geodesic curves. In [TT], the authors use a geometric model 
associated with Euler's elastica. The idea is to view the image as a surface 
and to minimize a linear combination of the curvature and arc length. More 
recently, the work in [S] , [TS] , and [21] provide state of the art results in the area 
of digital image inpainting. In [T3] and [^ a comprehensive survey of modern 
methods is given. 

We focus here on the model proposed by Bertalmi'o, Bertozi and Sapiro (here- 
inafter BBS), as presented in The model is based upon a transport equation 
analogous to the transport of vorticity in an incompressible fluid. The resulting 
model is a Navier-Stokes type time evolution. The approximated steady state 
solution of this model gives the corrected image. This model has several attrac- 
tive features. It is able to fill in regions surrounded by different backgrounds as 
well as regions that cross through boundaries, it is capable of handling arbitrary 
topologies, and as a third order PDE, it forces continuity conditions on both 
the image intensity as well as the gradient of the intensity across the boundary. 
The main drawback of the method is that it requires a large number of itera- 
tions (several thousands) to converge, making it not very competitive against 
noniterative fast- marching methods (e. g. [8]). 

In this paper, drawing motivation from the model presented in [5], we pro- 
pose to solve the third order PDE by using a different formalism, based on the 
minimization of the corresponding least-squares derived functional. The idea is 
to use for the minimization a steepest descent method with Sobolev gradients, a 
method that we found effective in reducing the computational time for different 
applications (e. g. [H]). The novelty and the challenge of the present approach 
is that we deal with a third order PDE, which is, in general, mathematically and 
computationally difficult to tackle. We establish in this contribution the math- 
ematical framework for the analysis of the method and show that the numerical 
performance is considerably improved compared to the original implementation 
of the model described in [5]. We thus offer a new possibility for using the 
Navier-Stokes model as an alternative to fast-marching methods, that remain 
unbeatable with regard to computing performance, but are still technical since 
they need a careful tuning of constants or functional parameters (see also [S])- 

In the method of calculus of variations, one often solves minimization prob- 
lems by computing the Euler-Lagrange equations which are used as a direction 
of descent for the energy. Explicit iterative schemes using the Euler-Lagrange 
equations require many iterations to converge due to a tight restriction on 



2 



the step size resulting from the CFL stabiUty condition. In our minimization 
scheme, we compute a gradient with respect to a Sobolev metric. In practice, 
this amounts to preconditioning the Euler-Lagrange equations for use in min- 
imization. The preconditioning in general allows for a larger time step when 
the resulting evolution equation is discretized in time. The Sobolev gradient 
method has proved to be very effective in a variety of image processing appli- 
cations, for example [10], [9], [22], and [25]. The focus of this work is to study 
the effect of preconditioning that is often associated with the Sobolev gradient, 
both in terms of the quality of the image and the computational efficiency. In 
contrast with the BBS model, we do not use a diffusion term, but instead ap- 
ply a smoothing operator to the Navier-Stokes gradient flow. An advantage is 
that our evolution equation is well-posed for all time and we obtain that the 
minimizer is across the boundary and can be obtained as an oj-limit of the 
evolution equation. The minimization of the energy functional with the Sobolev 
gradient falls into the more general category of image regularization [3], [3T], 
[55], [33], and [30] since we expect the interpolated image to be smooth due to 
continuity properties of the gradient. 

Although our results and applications are restricted to solving for the two- 
dimensional steady state arrived at by the Navier-Stokes model, our scheme is 
quite general. The ideas we present here can easily be adapted to computation- 
ally solve a vast variety of higher order partial differential equations. Since we 
treat the problem as a variational problem, our scheme can also be applied to 
existing variational models in image processing (which almost always use the 
Euler-Lagrange equations for minimization). An interesting, although not im- 
mediate, further application would be the use of the Sobolev gradient method 
for variational problems such as the ones presented in [2] . Finally, our scheme is 
cast in both an infinite dimensional Hilbert space setting and analogously in a 
finite dimensional finite-difference setting, making a strong connection between 
the theory of digital image analysis and applications. 

The outline of the paper is as follows. In section[21 we give a brief explanation 
of the BBS model as presented in [S] and make a connection to the Navier- 
Stokes model for incompressible fluids. In section [31 we describe the variational 
formulation from which we obtain our gradient flow in a Hilbert space setting. 
We show that the flow has a unique global solution and an w-limit at which 
minimization is achieved. In section |4] we obtain an expression for the Sobolev 
gradient. In section jS] we give the mathematical formulation in the discrete 
finite difference setting, and in section [6] we present the results from several 
examples which includes the resulting images as well as a comparison of CPU 
time, time step size, and the effectiveness of our scheme in solving the boundary 
value PDE. 

2 The Navier-Stokes model for image inpainting 

Here we give a summary of the Navier-Stokes model for image inpainting as 
presented in , also in [J , and [7| . A digital grayscale image can be thought of 
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as a set of m by n data points taking values between and 255. The resulting 
function / defined on the m by n grid is called the image intensity. An isophote 
of / is a level-line in this context. The smoothness of the image is represented 
by the Laplacian of the intensity. The direction of the level-lines of the intensity 
I at the (j, j)"* grid point is given by the direction that is perpendicular to the 
gradient of / at the point 

The idea of the Navier-Stokes model for image inpainting is to propagate 
the gradient of the smoothness of the images in the direction of the isophotes. 
Thus one wants to iterate the evolution equation 

I'{t) = V^/(<) • VA/(i) (1) 

or perhaps, by adding and anisotropic diffusion term, 

/'(<) = V-L/(t) • VA/(i) -f vV ■ {g\VI{t)\VI{t)) (2) 

to steady state solution 

V^/ • VA/ = 0. (3) 

In equation ([2|), a second term is added for its smoothing effects, where v > 
is small and g is a diffusivity function. 

The Navier-Stokes model for incompressible Newtonian fluids gives that the 
velocity field v and pressure p are coupled according to the equation 

v'{t) + v{t) ■ Vv(t) = -Vp{t) + iyAv{t) and V ■ v(t) = 0. (4) 

In two dimensions, the divergence free velocity field possesses a stream function 
ip so that V'^ip = V. The vorticity can be expressed a.s w = Ai/; and satisfies 
the advection diffusion equation 

w'{t) = -v{t) -Vwit) (5) 

when the viscosity is zero. The steady state must satisfy 

V^V • VAV' = 0. (6) 

This is the incompressible Euler equation for fluid flow. We refer the reader to 
for a background on the theory of the Navier-Stokes equation. 
Another way to say this is that the gradient of the stream function ip and 
the gradient of the Laplacian of the stream function must be parallel. Thus 
the analogy between the Navier-Stokes model and image inpainting is that the 
image intensity acts as the stream function in the Navier-Stokes model. In our 
approach, we consider an alternate to the method of Bertalmio as presented in 
[S] and minimize the norm of the left hand side using an evolution equation 
based on a Sobolev gradient method. 
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3 A variational formulation in the continuous 
setting 

In this section, we describe our scheme in a Hilbert space setting. The main idea 
is that we form an energy functional by taking the norm of the left hand side of 
equation ^ . We then search for a minimum of this energy using a gradient flow 
with a Sobolev gradient. We start by stating the definitions and results that 
we will need from the theory of Sobolev spaces. We then define the variational 
problem and study properties of the resulting functional. We argue that every 
critical point must be a solution of equation ([3|) . We define the Sobolev gradient 
in a general setting and, incorporating the boundary value into the gradient, we 
obtain a gradient fiow that has a unique global solution. We obtain a critical 
point of the energy as an w-limit of the gradient fiow. Since a critical point 
of the energy corresponds to a zero of the energy, we obtain that our energy 
functional does indeed posses a minimum. We discuss regularity properties of 
the minimizer. 

3.1 Sobolev spaces 

We review here the basic information regarding the theory of Sobolev spaces. 
The information provided here, as well as the notation, are taken from [1]. Let 
17 be a bounded open subset of . Recall that 



where for a multi-index a, D°'u denotes the a partial derivative of u taken in 
the distributional sense. Recall also that 



= the completion of {u e C""(r2) : \\u\\„r,p < oo} 
where C™(f2) is the space of m times continuously differentiable functions and 



with [| • Hp denoting the norm. is the closure of the the infinitely 

differentiable functions with compact support in il in W"^'P{fl). It is a result 
(Thm 5.37 of JL,), that under sufficient regularity conditions on the boundary 
of fl, the trace of D^u is zero for all u e and |a| < m. 

It is a result that for every open domain = W^™'P(r2). Since we 

wish to work in a Hilbert space setting, we Gx p ~ 2 and refer to the space 
H"^^^{n) as iJ™, = -ff^'^C^)' and L^{n) = L. 

Let Du = {D^w : 1 < |a| < k}. We can see that {(^^) : u G iJ™} is a 
closed subspace of L{il™') where £7™ — U\a\<m^a- Thus there exists a unique 
orthogonal projection P = P{m) from onto K^^^) : u G iJ™}. We men- 

tion this here as we will need this projection when constructing the Sobolev 



VF™'P(f7) EE {w g U\n) : D"u G LP{VL) for < |a| < m} 
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gradient in section We also note that is a closed subspace of iJ"*. Thus 
there exists a unique orthogonal projection of i?™ onto H™ which we denote 
hyPo^Poim). 

Suppose that fl is bounded and has C™ boundary dfl. It is a result, see [T], 
that for u G i?™, the boundary traces of u can be defined as fa{u) = Dau\gQ 
for < |a| < m. Further is in L^(dn). If u e i/^", then = for ah 

a. Finally, we remind the reader the part of the Rellich-Kondrachov theorem 
that we will need. 

Theorem 1. Suppose that 51 C is a bounded domain with a smooth boundary. 
The embeddings 

and 

are compact when to > 2. 

Here C^(f2) denotes the j times continuously differentiable functions on fj. 



3.2 Defining the minimization problem 

Assume that is a bounded open set in and has smooth boundary. For 
u G , define F so that 

F{Du) = V^u ■ VAu. (7) 

Now if u e , then D°'u £ when |a| — 1. Since is compactly embedded 
in C(f2), for u G H^, Vm is a bounded continuous function, and |Vu| < c||u||//3 
for some constant c independent of u. | • | denotes the sup norm here and in the 
rest of the paper, || • || denotes the norm, and || • ||_f/m the i?™ norm. Thus 

\\F{Du)\\ = (^j^ |V^u • VAu|2^ ' < \Vu\ ||VAu|| < c \\u\\]j, (8) 
and we see that for u G ^ F{Du) G and thus the energy 

Eiu) . (9) 

is well defined. To reduce notation, we write uq = P{)U, where Pq is the orthog- 
onal projection of onto Hp, . 

In the model of BBS, u G is sought so that F{Duo) — and u = g on 
the boundary of fl for some predefined function g. We assume that g is on 
the boundary of fl. 

First we show that E is differentiable then explain why a minimizer of 
E would correspond to a zero of F o D. 

Theorem 2. E : — )• M as defined in equation ([9]) is Frechet differentiable. 
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We compute here the first and second Frechet derivatives of E. We regard 
the first Frechet derivative of at m £ as a bounded hnear operator on 
and the second Frechet derivative of i? at u G as a bounded and symmetric 
biUnear operator on H^. 

Proof. Wc compute a Frechet derivative of F as given in equation (O to obtian 

F'{Du)Dh = \7^h-VAu + V^u-VAh. (10) 

One can also compute F" {Du){Dh, Dk) in order to see that F o D is 
Frechet difFerentiable on H^. Next we can compute the first and second Frechet 
derivatives of E to see that 

E'iu)h = {F'{Duo)Dho, FiDuo)) (11) 

and 

E"{u){h,k) = {F'iDuo)Dho,F'{Duo)Dko) + {F"{Duo){Dho, Dko), F{Duo)). 

E being Frechet difFerentiable follows from the differentiability of F. 

Notice from equation ([TU]) that F'{Du)Du = 2F{Du). From equation pT|) . 
we see that 

E'{u)uo = {F'{Duo)Duo,F{Duo)) = 2\\F{DuoW = ^^Eiu). (12) 

□ 

Notice also from equation (fT2|) that E'{u)uq — 2\\F{Duo)\\'^ . Hence, in order 
to solve the boundary value problem, it suffices to find u e H'^ so that u = g 
on on and E'{u)h = for all h G Cf^{n). 

3.3 The Sobolev gradient 

The idea of Sobolev gradients [32], is that given a function (p that is every- 
where defined on a Hilbert space H, we can represent the Frechet derivative of 
E aX u using a member of the Hilbert space. This is due to the Riesz represen- 
tation theorem as if is C^, then (t>'{u) is a bounded linear functional on H. 
Thus there exists a unique element of i?, which we denote by Vi?(u), so that 

(t)'{u)h = (/i, V0(u))ff for all h(^H. (13) 

One then considers the gradient flow 

z(0) = xoe H and z'{t) = -V(/)(z(i))- (14) 

It follows that z defines a path along which the energy is decreasing. Assuming 
one has existence of such a z, then it is important to study the asymptotic limit 
as it is in this limit that we hope to achieve a critical point. 

For this problem we consider a gradient with respect to the H'^ inner product. 
In order to obtain an w-limit however, we need to map this gradient into a higher 
order space. In this section, we give two results from [32] that wc will modify 
for this work and define our evolution equation. 



7 



Theorem 3. Suppose (f> is a nonnegative valued function on a Hilbert space 
H . Define the gradient of (j) at u as in equation (jl3p . IfV(j) is a locally Lipschitz 
function from H to H, then for xq ^ H the gradient flow pi]) has a unique 
solution for all t > 0. 

Theorem 4. Suppose, under the conditions of theorem that z is a global 
solution of (jl4p with energy (p. Further suppose that (f>'{u)u > for all u G H , 
then the range of z is a hounded subset of H . 

The main idea of this theorem is that 
(Mint) = {z{t),z'{t))H - -{z{t),VH^{zmH - -cl^'{z{t))z{t) < 

by assumption. This imphes that t \\z{t)\\'jj is non increasing and thus the 
range of z is bounded in H. 

Since E, as defined in equation ([9]), is Frechet differentiable, E'{u) is a 
bounded hnear functional on . Thus for each u G , there exists a unique 
member of so that 

E'{u)h = {h,VH'iE{u))H^. 
For X e and fc > 3, consider the mapping from to M given by 

y -> {y,x)Hi- 

One can see that this is a bounded linear mapping. Thus there exists a 
unique element Mk^ x G Hq so that 

{y,x)m = {y,Mk,o x)Hk 

for all y G Hq. Now consider the operator Mk^o- In [T7], we obtained that this 
operator is bounded from Hq to Hq with norm less than or equal to one. 
We define the gradient flow 

z(0) = uo and z'{t) = -MkflPQ^ mE{z{t)). (15) 

where uq G H'^ and uq = g on dfl. Here Pq is the projection of onto Hq. 

3.4 Global existence, uniqueness, and asymptotic conver- 
gence 

In this section we give results regarding the global existence and uniqueness of 
the gradient flow ((T5t . We show that a unique solution exists for this system and 
that the range of z is bounded in H'^. This allows us to extract a subsequence 
from the range of z which converges to a zero of the gradient thus giving the 
existence of a stationary point and a solution to the minimization problem. We 
follow the developments in chapter 4 of ^32 closely. 

Theorem 5. The gradient flow given in equation (jlSp has a unique global so- 
lution z e cH[0,oo),iJ'^'). 
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We follow the proof of theorem [3] presented in [35] with only a slight modi- 
fication to account for the presence of MkflPa in our gradient flow. 

Proof. In order to show that the gradient flow given in (|15p has a local solution, 
we need to show that the gradient is locally Lipschtiz from H'' to H'^. We 
showed in theorem [5] that E is on . Note that for any function </) 
defined on a Hilbert space H, with gradient as in ([T3| . there exists a constant 
Cu and a ball B about u so that if u g _B then 

V0(m) - V(/.(w))h| = |(<^'(u) - < Cu\u ~ v\H\h\H. 

Thus if we take h = \74'{u) — V0(t;), then we see that 

\V(f){u) - V(f){v)\H <Cu\u- v\h. 

Hence if is a function defined on a Hilbert space, then the gradient of (j) 
as defined in ((T3)) is locally Lipschitz. From this we conclude that V3E : — > 
is locally Lipschitz. Mk.oPoV hsE is locally Lipschitz from H'^ to as for 

u, V e H^, 

\\MkfiPo{V H^E{u) -V h-.E{v))\\h^ < \\Po{VH.Eiu)-VH.Eiv))\\H3 < 
\\{\7h3E{u) - \7h3E{u))\\h3 < Cu\\u - v\\jf3 < c„||w - vWh^. 



We conclude that the system ([TSj) has a local solution. 

Let T be the largest number so that a solution exists on [0, T). In order to 
show that the flow has a global solution, it suffices to show that limf^y- z{t) 
exists. This holds true provided that there is a constant m so that 

\\z{a) ~ z{b)\\ffk <m (16) 

for all < a, & < T. Recalling that Pq is symmetric with respect to {■,-)h^ and 
z'{t) e Hq and hence in the range of Pq for all t, we have that 

' , V 

\\z'\\h>' < (using Cauchy-Schwarz) 
(6-a) / \\z'\\l,={b-a) I {z',z')h^ =-{b-a) [ {z' , MkflPo^ h^E{z)) < 

z',Po\/H3E{z{t)))H3^~T f {Poz\\/H3E{zmH3 = 

J a 

T I {z',\7H^E{z{t)))H3 ^-T I E'{z{t))z'{t)^~T I {E o z)' < TE{z{0)). 



b rb rb 

l\\2 




Thus equation ([T5)) holds true and we have that the flow ([T5|) has a global 
solution. Uniqueness follows from the fact that the projected gradient is a locally 
Lipschitz function and the fundamental theory for existence and uniqueness of 
ODE'S. □ 
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We show next that the flow (|15l) has an w-limit which is a minimizer of the 
energy and a solution to the boundary value problem. 

Lemma 1. There exists an unbounded sequence of numbers {tn}n>i so that 
Mk,oPi^V ukE{z{tny) converges to zero in . 

Proof. Using the same analysis as in the previous proof, we note that 

/■oo /"OO 

/ UhflP^V H-<^E{z)\\]j, <- {Eoz)' < E{zm. 
Jo Jo 

From this it follows that 

poo 

\\Mk.oPoVH-E{z)\\]j, 



is bounded and hence the conclusion follows. □ 

Lemma 2. The range of z is bounded as a subset of . 

The proof is similar to the proof of theorem |4] given in 132; with only a slight 
modification. 

Proof. Let h{t) = .5\\z(t)\\'j^k and note that 

h'{t) = {z{t),z'{t)}H>^ = -{zit),Mk,oPo^H'^Eiz{t)))H^ = 
-{Pozit),VH^E{z{t)))H^ = -E'{z{t))P^z{t) - -AE{z{t)) 

the last equality following from (|12p . Thus h' is never positive and hence h is 
bounded. □ 



Theorem 6. Suppose k > 5 in equation (|15p . Then there exists a sequence of 
unbounded numbers {tn}n>i so that {z{tn)}n>i converges in to u ^ C^(f2) 
and E'{u)h = for all h e C^{Vt). 

Proof. Since the range of z is bounded in H'^ for fc > 3, using the compact 
embedding of into , there exists a sequence {t„}„>i so that {z{tn)}n>i 
converges to u G . 

Using theorem[Tl we can assume {t„}„>i is chosen so that {Mk.oPo"^ H'^E{z{tn))}n>i 
converges to zero in . Using the regularity properties derived in theorem [5l 
Mk,oPo'^H^E{z[tn)) converges to in . For h € C^, 

E'{u)h^ {h,VH3E{u))H^ = 

{Poh,VH^E{u)}HS = {h,PoVH^E{u))H'. = 

{h,Mk,oPoVH-,E{u))H^ =0. 
Since is embedded in C'^iQ), u e C\n). □ 
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We also note that z{t) agrees with g on the boundary of ft for all t. Recall 
that z has range in H'' which is compactly embedded in C^(f2). Let f{t) = 
zit)\on- Then / is and l/'(t)| = \Mk,oPoVH3E{z{t))\9n\ = 0. Since 

\f{t)~9\ = \fit)~fm< t\f'\ = 0, 

Jo 

the assertion follows. The norm | • | here indicates the sup norm. 

Regarding the uniqueness of the solution of the boundary value problem 

V^u ■ VAu = in and u = g on dQ, 

we refer the reader to [5j where an example is given to illustrate that the solution 
of the boundary value problem is not unique. It is also argued that some level 
of non uniqueness is required for the inpainting problem as one would like to be 
able to 'choose' the best interpolation for the inpainting problem. 

4 An expression for the gradient 

Here we give an expression for the Sobolev gradient. We first derive the Euler- 
Lagrange equation for our energy functional. We observe that the Euler-Lagrange 
equation can be viewed as a gradient of the energy obtained with respect to an 
inner product. We then obtain an expression for the Sobolev gradient us- 
ing the projection P discussed in section [3. II and compare the Sobolev gradient 
with the Euler-Lagrange equation making a key distinction between the two by 
comparing the regularity properties of the resulting gradients. 

4.1 The Euler-Lagrange equation for E 

Recall the expression for E'{u) as given by equation PT|) . In order to obtain the 
Euler-Lagrange equation, we need first to find some vector v whose components 
are in so that 

{F'{Du)Dh,F{Du)) ^ {Dh,v). (17) 

For this purpose, let y ^ and recall the expression for F' {Du)Dh as given 
by equation pH)) to obtain 

{F'{Du)Dh, y) = (V-L/i • VAw + V^u ■ VAh, y). 

Since Sl^h ■ VAu ~ — V/i • V^Au, we have that 

{F'{Du)Dh, y) = (V/i, -y • V-^Au) + (VAh, y ■ V^u). 

Here, when c is a scalar, c • (^ = . 

From this calculation we have that the nonzero components of v are Va-^ = 
-F{Du) ■ V-^Au and = F{Du) ■ u so that 

{Dh,F{Du)) = (V/t,w„,) + {VAh,v,,,). 
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With this notation, under further regularity assumptions on u, the Euler-Lagrange 
equation for E becomes 

VELEiu) = D*v (18) 

where D* denotes the adjoint of I? as a closed and densely defined operator. 
Note that for all h e L'^ 

E'{u)h = {h,VELE{u)). 

Thus the Euler-Lagrange equation, when defined, can be viewed as an gra- 
dient for E. 

4.2 An orthogonal projection 

Recall from section 13.11 that P is the orthogonal projection of L^(r2'^) onto 
{(m) ■ ^ ^ with Dh denoting all partial derivatives of h up to order k. 
Recall also that the expression for E'(u) as given by equations ^TU\\ and pT|). 
Then we have the following 

E'{u)h = {F'{Du)Dh,F{Du)) = 




where n(p = x. An expression for P, given in [32 , is 

f {I + D*D)-^ D*{I + DD*)-^\ 
-\D{I + D*D)-^ I-{I + DD*)-^) 

where D* denotes the adjoint of D when viewed as a closed and densely defined 
linear operator on i^(f2). From this formula we see that 

Vni-Eiu) = p(^_^ = D*{I + DD*)-'^v. (19) 

We note that if u satisfies additional regularity so that v is in the domain of 
D*, then 

D*{I + DD*)-^v^ {I + D*D)-^D*v. (20) 

Thus the Sobolev gradient can be viewed as a smoother gradient than the Euler- 
Lagrange equation which is the descent direction used in most optimization 
problems. In case k = 1, {I + D* D)~^ is the resolvent of the Neumann Laplacian 
(/ — A)~^. For fc > 1, a nice argument is given in [S] to show that (/ -I- D*D)~^ 
is equivalent to (/ — A)^'^. Thus we see that using a gradient in a higher order 
Sobolev space is equivalent to preconditioning the Euler-Lagrange equations 
using the resolvent of the Neumann Laplacian as the preconditioner. A detailed 
discussion on the preconditioning resulting from the use of gradients in high 
order Sobolev spaces is presented in the Appendix [Al 
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5 Finite dimensional implementation 



We discretize, using finite differences, equation ^T9\\ or equivalently ((20)l since 
in a finite dimension space tfie two expressions are equal. Here u denotes the 
image intensity, a real valued function defined on a rectangular domain i? of ffi.^ . 

5.1 The inpainting domain 

The inpainting domain is specified by the user. This can be done by marking the 
regions to be inpainted using a photo editing software such as Gimp or Paint. 
We then determine the inpainting region, fi, and the inpainting region plus its 
boundary, il'. We consider u{i,j) to be in ft' if u{k,l) G fl for some {k,l) with 
< |i — A:| + |j — Z| < 2. We label uq to be u restricted to il and u' to be u 
restricted to fl'. uq and u' are converted into vectors using the natural ordering 
+ > (ij) then (ij + l) > 

5.2 The finite difference operators Di, D2, and A 

Let Di, D2, and A be matrices defined on the interior of R so that 



for all G fl- We use second order centered differences to approximate 

the first partial derivatives and the five point discretization of the Laplacian. 
Define also Di A = Di*A and I?2 A — D2 * A. In order to enforce the boundary 
conditions, we eliminate the rows of each of these operators corresponding to 
pairs G i? — and columns corresponding to pairs G i? — il'. Finally 
we define (/ — A)^^ : — > to be the resolvent of the Laplacian. 

5.3 Discretizing the gradient 

Here we discretize the Sobolev gradient. Note that in the finite dimensional 
case, equation pp]) holds true. First we discretize F{Du) as 



where for two vectors x and y, x y denotes the vector consisting of their ele- 
ment wise product. Next we discretize the Euler-Lagrange equation as follows. 
Discretize v given in ([TBI) so that 



D2Uij = ^(uij+i - Ui,j-i), and 



F{Du) = D2U' * DiAu' - Diu' * DaAu' 



(21) 
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Then the discretized Euler-Lagrange equation for this problem becomes 

QEL = D[v^, (1) + D'^v„, (2) + ^{D[v^, (1) + D'^v^, (2)) (22) 
where D[ denotes the adjoint of Di. 

5.4 Smoothing the Euler-Lagrange equation 

In order to use the Sobolev gradient as given in equation (1^ . we need to 
discretize the operator {I+D*D). Once we have a discretization of this operator, 
we obtain the discrete Sobolev gradient by solving the equation 

(/ + D*D)gsob = gsL 

Since D involves partial derivatives or order up to three, a discretization of 
this operator would force us to solve a very complex linear system which can 
be computationally very time consuming. A better approach is to obtain an 
equivalent operator to use for smoothing the Euler-Lagrange equation. A good 
choice turns out to be (/ — A)'^. A good explanation for this is given in [5]. 
Thus in order to solve for the Sobolev gradient, we solve the linear system 

{I~Afgsob^gEL. (23) 

This is computationally the most expensive part of the algorithm. To solve the 
linear system, a Cholesky factorization of (/ — A) is performed to obtain an 
upper and lower triangular matrix. Then the linear system 

{I-A)x = y. 

is solved three times using the factorization in order to apply the smoothing 
operator three times. This is convenient method since homogeneous Dirichlet 
boundary conditions naturally apply at each step of the factorization. 

5.5 Discretization in time 

In order to discretize the system (jlSp in time, we use an explicit Euler scheme 
with a locally minimizing time step. Thus we generate a sequence of images 
according to 

uo{n + 1) = uain) - tngsob{n). 

We remark that several options are available for time discretization. One can 
choose between an explicit or implicit scheme. One can also use a fixed time step. 
For our purpose, we chose the explicit scheme since computing the gradient here 
is more expensive than computing the value of the functional. Thus an implicit 
scheme that would require us to compute the gradient several times at each step 
is a lot more expensive than an explicit scheme with a line search algorithm to 
compute locally the optimum time step. 
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6 Results 



We present here the results from our simulations. The implementation of the 
discrete operators derived in the previous section using finite differences is quite 
straightforward if a language with vector programming capabilities is used. The 
algorithm is programmed in Matlab and all operations are written in matrix 
form. We run our simulations on an IBM laptop with 1.5 GB RAM and 1.4GHz 
Pentium M processor. In each case the starting point for the simulation was the 
solution of Laplace's equation 



which was solved by using successive over relaxation method. The boundary 
condition function g takes values from the neighboring non-inpainted domain. 
In the examples we present, we normalize all images so that the maximum pixel 
value is one. This choice for the initial guess of the solution looks natural if 
the Navier-Stokes analogy is used: equation ((M)) then represents an irrotational 
fluid field, i. e. with zero vorticity oj = A?/; = 0, that will be iterated to find 
a nonzero vorticity (source) term that is present in the model. In the method 
of BBS, the initial condition is set after performing several steps of anisotropic 
diffusion following ([2]). 

For comparison purposes, we also implemented the method of BBS as is 
described in [4j . We incorporate slope limiters for convection term to propagate 
fronts with sharp variation and apply every 50 time steps a few iterations of the 
Perona-Malik scheme for anisotropic diffusion to avoid spurious oscillations: 



As many of the specific details such as the term used for anisotropic diffusion 
are not available in the original paper of BBS, an exact reproduction of their 
results proved to be difficult. As a consequence, the comparison presented in the 
following should be regarded with caution. Nevertheless, reported convergence 
properties of the method of BBS are recovered by our implementation. 

6.1 Comparison with the original implementation of the 
Navier-Stokes model by Bertalmio, Bertozi and Sapiro 

We perform a full quantitative comparison displaying convergence curves using 
our minimization procedure and the Navier-Stokes model of BBS. In the orig- 
inal model it is proposed to solve the equation V^u • VAu = by evolving 
this PDE in time adding anisotropic diffusion. Typical methods coming from 
Fluid dynamics literature are used, and therefore, several delicate and technical 
choices are necessary: the choice of the limiters for the discretization of convec- 
tion terms, the relaxed Laplacian solvers for the vorticity needing the tuning 
of a constant, the choice of how often and under which form to introduce the 
diffusion term to avoid spurious oscillations. Besides, no convergence analysis 
is provided to state that a solution of the PDE is obtained within a reasonable 



Am = and = g 



(24) 



I'{t) = V ■ c{x, y, t)VI{t) where c{x, y, t) = e^H'^-f "/^ 



(25) 
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number of iterations. A rough estimation of the necessary steps for convergence 
is given in [8], using the CFL time-step restriction for the exphcit scheme used 
by BBS. The number of steps to reach a stationary state is estimated as large 

1/2 

as N^-^^i^, which results in several thousands of iterations for typical test cases. 
Our implementation of the original method displays convergence properties of 
this order of magnitude. 

Our method does not display these shortcomings related to the practical 
implementation: there are not constants to be tuned and no advanced knowl- 
edge of numerical methods for Navier-Stokes equations is required. Since we 
use a minimization technique, we can choose a locally minimizing time step 
to ensure that the energy is minimized. In addition, by preconditioning the 
Euler-Lagrange equations, we remove the need to apply anisotropic diffusion. 
In contrast with the method of BBS, the results for inpainting with and without 
anisotropic diffusion were nearly identical for the test cases we tried. 

We first evaluate convergence properties of our method, compared to the 
method of BBS, for the test case displayed in Fig. [2] Two minimization schemes 
were used, based upon the Euler-Lagrange equations and the Sobolev gradient 
using just one application of the preconditioner (also refered to as the gradi- 
ent), respectively. For the purpose of demonstrating the effect of preconditioning 
on image inpainting, we found it sufficient and reasonable to apply the precon- 
ditioner (J — A)~^ just once. Thus in the later images in order to save CPU 
time, we used just one application of preconditioning. 

We plot in Fig. [T] the reduction in ||V-'-it • VAu|p using the method of 
BBS and our energy minimization method. We can observe that our method is 
quantitatively closer to the solution of the PDE. With our minimization scheme 
we were able to achieve a decrease of two or three orders of magnitude for most 
cases. As our method is a gradient descent method, we expect a linear rate of 
convergence. However, we do note that by using the Sobolev gradient in place of 
the Euler-Lagrange equations, the convergence to a zero is slightly accelerated. 

In the same figure [1] we plot the "error" e = max|u„+i — where Un 
represents the current image. This quantity is used as stoping criterion for the 
iterative process since it seems to be the best measure for convergence in image 
processing, as it is generally a valid measure of convergence in numerical anal- 
ysis. We should note in passing that, compared to Fluid dynamics applications 
where e < 10~^ is required for convergence, for image processing applications 
no apparent change in the image is observed starting with e < 10"'*. 

An important question related to convergence is the condition number of the 
algorithm. Although the condition number of the least squares problem when 
one considers minimization with the Lagrange equations is relatively high, we 
demonstrate, that this number is significantly reduced when we precondition 
the problem and take the gradient in a higher order Sobolev space. Since these 
results concern more general numerical analysis considerations, we report them 
in the Appendix \X\ 

The images in Fig. [2] were ran until an error tolerance of 10~^ was achieved. 
Using the method of BBS, this took 9282 iterations in a CPU time of 67.28 
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Figure 1: Left: Comparison of the norm of the residual (energy) for the different 
methods. Right: The error as defined by max |u„+i — Un\ for the different 
methods. N-S denotes the Navier-Stokes method following the implementation 
of Bertalmfo et al. E-L denotes the iterative method using the Euler-Lagrange 
equations. denotes the Sobolev gradient. 



seconds. Using energy minimization with the Euler-Lagrange equations, this 
took 1014 steps in a CPU time of 29 seconds. Using one application of the 
preconditioner, convergence was reached in 742 iterations in a CPU time of 32 
seconds. Using three application, convergence was reach in 942 step in about 
40 seconds. We remark that by preconditioning we were able to reach the error 
tolerance using a lesser number of iterations. The CPU time was comparable 
for the case using the Euler-Lagrange equations and the Sobolev gradient gra- 
dient, as solving a linear system at each time step adds to the processing time. 
However, one could take further measures to decrease this additional time by 
implementing a different algorithm for solving the linear system. 

Qualitatively the effect of preconditioning is apparent in the bottom two 
images in figure[2j While the image resulting from the Euler-Lagrange equations 
has some noise where the inpainting was performed, the image resulting from 
the Sobolev gradient is smooth in the inpainting region. Also note that the 
edge between the face and the background is correctly placed. In this sense, we 
remove the need to apply anisotropic diffusion. 

In this image the result we present for the BBS method is qualitatively 
worst than the result presented in the original paper of Bertalmfo et al. [4] . We 
attribute this to the fiuctuating error as presented in figure [T] the result for this 
problem initially improved, however the subsequent images did not correctly 
inpaint the eye. This behavior is certainly related to the careful choice of the 
parameters in the original implementation (time step, when and how to apply 
the anisotropic diffusion) . Since the error analysis is not available for the method 
of BBS, we didn't attempt to find a deterministic choice for these parameters. 

Since the full image of girls, first appeared in [4], is currently used to test 
inpainting methods, we present in Fig. |3] our results for this test case. The 
inpainting region is correctly filled, with a better image quality obtained when 
the Sobolev gradient method is used. 
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Figure 2: The resulting images for the simulations. Top left: the defective 
image, Top right: Navier-Stokes inpainting, Bottom left: Energy minimization 
with the Euler-Lagrange equations, Bottom right: Energy minimization with 
the gradient. 

6.2 Qualitative comparison with other methods 

In an effort to compare our results to other methods used in image inpainting, 
we use here the image |4] taken from [8]. We performed our inpainting algorithm 
just on the gray scale image as the issues of edge recovery and smoothness of 
the interpolation can be determined from this image. 

In comparison with the method of Bornemann and Marz [5j , the construction 
of the iris is comparable. In their result a better quality image is achieved as 
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Figure 3: Inpainting of the girls. Compare with the original papers f4| and [8]. 
Top left: original image, top right: the solution of Laplace equation, bottom left: 
the solution using Navier-Stokes inpainting, bottom right: energy minimization 
with the gradient. 



the texture of the eyebrow is interpolated. Also the CPU for this example was 
about 60 seconds. In [8 a reported CPU time of 0.4 seconds was reported. It 
goes without saying that our iterative method could not compete in terms of 
efficiency with a non- iterative fast-marching method, such that proposed in [S]. 
Nevertheless, the method of Bornemann and Martz needs a careful tuning of 
four parameters, which is not the case for our method (see also [6]). 



6.3 Image interpolation 

In this section, we apply our algorithm to perform image interpolation. We 
start with the two images given in figure [5] and interpolate the bird to 160 by 
152 pixels and Lena to 200 by 256 pixels using nearest neighbor interpolation in 
figure [SI The interpolation is carried out as follows. The image is first expanded 
to the dimensions of the larger image where the newly added pixel designate 
the inpainting region. As with inpainting, we first solve Laplace's equation then 
apply the inpainting algorithm to further sharpen the image. 

For these two images the effect of using anisotropic diffusion in the method of 



19 




Figure 4: Inpainting of the eye as was done in [5|. Top left: original image, top 
right: the solution to Laplace's equation, bottom left: the solution using Navier- 
Stokes inpainting, bottom right: energy minimization with the gradient. 




Figure 5: Image interpolation: original image enlarged using nearest neighbor 
interpolation. Left: bird original size 80 by 76 pixels. Right: Lena original image 
50 by 64 pixels. 

BBS is evident. The smoothing effect of anisotropic diffusion does not allow one 
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to recover the detail of these two images. In the energy minimization method, 
since we avoid using anisotropic diffusion, we can recover sharp edges as well as 
some of the finer detail. For example in the image of the bird, notice that in 
the energy minimization method using both the Euler-Lagrange equation and 
the Sobolev gradient, the eye of the bird was actually recovered. The effect 
of preconditioning here is evident in that it minimizes some of the artifacts 
associated with image interpolation. For example comparing the beak of the 
bird in the image resulting from cubic interpolation and minimization with the 
Lagrange equation vs minimization with the Sobolev gradient, we see that the 
latter results in sharper edges. Similar effects are evident when we examine the 
interpolation of Lena. In particular, the smoothing and sharper edges are most 
visible in the upper part of the hat of Lena. 




Figure 6: Resulting images for image interpolation. Interpolation to twice as 
many pixels. Top left result from the method of Bertalmi'o et al. , top right result 
from bicubic interpolation, bottom right result from minimization using the 
Euler equations, bottom right result from minimization using the gradient. 
Notice artifacts are not present when using the gradient (see beak of bird 
for example). 

We remark here that in practice the assumption that u be continuous 
on the entire domain Q is too strong. In image inpainting an edge represents 
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Figiirc 7: Resulting images for image interpolation. Interpolation to four times 
as many pixels. Top left result from the method of Bertalmio et al. , top 
right result from bicubic interpolation, bottom right result from minimization 
using the Euler equations, bottom right result from minimization using the 
gradient. 

a discontinuity, thus in theory this important case would be excluded from 
our analysis. The idea of using a gradient in a higher order Sobolev space is 
to increase the smoothness of the interpolation, and this is supported by the 
theory as when one computes a gradient in a higher order Sobolev space, we see 
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that the minimizer should also satisfy further regularity conditions. In several 
of our cases, using a Sobolev gradient reduced artifacts associated with image 
inpainting and image interpolation while preserving edges (see images [51 El and 
[Hand the associated discussions. However, some care should be taken in practice 
as over smoothing would result in blurred reconstruction of edges. This artifact 
appeared in some or our results (see for example H]). In principle, the analysis 
of this problem could be studied and carried out in a lower order Sobolev space 
or in the space of functions of bounded variation which is the traditional space 
used for an analytical study of problems in image processing. We expect that 
the estimates required for such an analysis would be of a similar nature to 
the estimates required for existence and uniqueness results of the Navier-Stokes 
equations and hence beyond the scope of this paper. 

7 Summary and final discussion 

In this work we derived a method to solve the digital image inpainting model 
that is based on the model of the Navier-Stokes equations for an incompressible 
flow. We used a variational formulation and a gradient flow to converge to the 
stationary solution. The flow is based on a gradient coming from a Sobolev 
norm, instead of the generally used norm, which in practice is associated 
with preconditioning the Euler-Lagrange equations. We found that through the 
use of the preconditioner, we can favorably alter properties of the interpolation 
by reducing noise and offering a more accurate interpolation of the edges. 

Compared to the original Navier-Stokes model proposed by Bertalmio, Bertozi 
and Sapiro (BBS) 5 , our algorithm displays better convergence properties, with 
similar qualitative image inpainting results. From the practical point of view 
of computer implementation, we have eliminated technical aspects related to 
the use of numerical methods for Navier-Stokes equations (slope limiters, dy- 
namic relaxation for Poisson equation, anisotropic diffusion steps, etc). We also 
mention that the new algorithm does not require to tune constants, as in non- 
iterative methods [5] . Using a straightforward finite difference implementation, 
we demonstrated, through various examples for image inpainting and image 
interpolation, that the algorithm is effective for such problems. 

We finally should like to emphasize that the mathematical theoretical frame- 
work that we developped in this paper is of more general interest and could be 
applied to other variational formulations used in image inpainting, such as the 
one presented in [TH]. Such new applications of the Sobolev gradient appear as 
a natural extension of the present work. 
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A Condition number 



We give here two precise definitions for the condition number and show that by 
Sobolev preconditioning, the relative condition number is, as shown in figure [H 
improved. 

Following the definition given in ^S] , the following definitions for the condi- 
tion number of a function / from an arbitrary normed space X to an arbitrary 
normcd space Y is used. 

Definition 1. The relative condition number of a differentiate function f : 
X ^ Y is defined as 

\\f{x)\\L(X.Y)Mx 
\\fi^)\\Y 

Definition 2. The absolute condition number of a function f : X ~^Y is given 
by 

y \\fix + h)~fix)\\ 
.= hm ..P||,|,,<, ^ . 

||/'(a;) is the norm of f'{x) when viewed as a bounded linear operator 

from X to y for a differentiable function /. In case X = M" is equipped with 
the Euclidean metric and y = R, then f'{x) is just the Jacobian of / at x. 

Note that if y = R, then f'{x)h = (h, Vxf{x))x for each x ^ X, following 
the definition of the Sobolev gradient. Thus ||/'(a;)||L(x r) — sup{f'{x)h : h G 
X, \\h\\x = l} = \\Vxf{x)\\x. 

In the finite dimensional setting, iJ*^ is the space of all n-dimensional vectors 
equipped with the metric (u, v) [{k — (u, (/ — A)'^w)Rn . Here A is a discretization 
of the Laplacian. Let / : H'^ — ^ R, then using the definition of the gradient 

f{x)h = {h, V/(a;))R. = {h, {I - A)'^Vff./(x))K. Vh £ R" . 

Thus 

VH.f{x) = {I-A)-'^Vfix) (26) 

and 

I|V^./(x)||h. = ((/-A)'=(/-A)-^V/(a;),(/-A)-'=V/(a:)) = (V/(a:), (/-A)-^V/(a;)). 

We make two remarks. First, we see that by equipping the domain of the 
problem with the H'^ norm in place of the Euclidean norm, we precondition the 
problem so that the absolute condition number of the problem as defined above 
is reduced by a factor of |(/ — A)^^!''. 

This effect could be easily illustrated if the Fourier representation of the 
gradient is used. If the ordinary gradient is sampled on a rectangular, regular 
grid as: 

V/(x,) = ^5e*^''-^ (27) 
p 
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then, following (P^ . the gradient will have a similar representation (same 
basis functions) with modified Fourier coefficients: 



Since /c > 1 we infer that ||V/ffc/(a;)|p — ^ is reduced and therefore the 

condition number is improved. The practical consequence is that diminishing 
the corresponding wave numbers in the Fourier decomposition results is attenu- 
ating the most oscillating modes of the solution. This could explain the smooth- 
ing effect of Sobolev gradient methods (see also 16 ). Second, using the above 
reasoning, if / is real valued and X = H'^, then \\f'{x)\\L(^x.Y) — ||Vx/(a;)||x 
and hence the relative condition number is given by 

{Vf{x)Al^A)-'^Vf{x))\\x\U 
1/(^)1 

In the examples we present, we normalize all images so that the maximum 
pixel value is one and thus use the relative condition number as a measure of 
how well conditioned the problem is. 

In figure [51 we give a plot of the relative condition number for the case 
corresponding to image [2j We compare the relative condition number of the 
problem with and without preconditioning. The relative condition number is 
lowest for H'^, achieving a reduction of a factor of 4 compared with the Lagrange 
equations. 

Another possibility to deal with the condition number issues is to refor- 
mulate this problem as a first order system. This method is currently under 
investigation in |23j and could be possibly applied to the inpainting problem. 
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